\chapter{Iterative methods}
\label{sec:iterative_methods}

The following sections deal with the methods of Iterative Boltzmann Inversion (\ibi) and Inverse Monte Carlo (\imc).

In general, \ibi and \imc are both implemented within the same framework. Therefore, most settings and parameters of those methods are similar and thus described in a general section (see sec. \ref{sec:iterative_methods_imc}). Further information on iterative methods follows in the next chapters, in particular on the \ibi and \imc methods.

\begin{figure}[h]
   \centering
   \includegraphics[width=7cm]{usage/fig/flow_ibi.eps}
   \caption{\label{fig:flow_ibi}Flowchart to perform iterative Boltzmann inversion.}
\end{figure}

\section{Iterative workflow control}
\label{sec:iterative_workflow}

\begin{figure}[t]
   \centering
   \includegraphics[width=7cm]{functionality/fig/flowchart.eps}
   \caption{\label{fig:flowchart}Block-scheme of the workflow control for the iterative methods. The most time-consuming parts are marked in red.}
\end{figure}

Iterative workflow control is essential for the \ibi and \imc methods.

The general idea of iterative workflow is sketched in fig.~\ref{fig:flowchart}. During the global initialization the initial guess for the coarse-grained potential is calculated from the reference function or converted from a given potential guess into the internal format. The actual iterative step starts with an iteration initialization. It searches for possible checkpoints and copies and converts files from the previous step and the base directory. Then, the simulation run is prepared by converting potentials into the format required by the external sampling program and the actual sampling is performed. 

After sampling the phasespace, the potential update is calculated. Often, the update requires postprocessing, such as smoothing, interpolation, extrapolation or fitting to an analytical form.

Finally, the new potential is determined and postprocessed. If the iterative process continues, the next iterative step will start to initialize.

\subsubsection*{How to start:}

The first thing to do is generate reference distribution functions. These might come from experiments or from atomistic simulations. To get reasonable results out of the iterative process, the reference distributions should be of good quality (little noise, etc).

\votca can create initial guesses for the coarse-grained potentials by boltzmann inverting the distribution function. If a custom initial guess for an interaction shall be used instead, the table can be provided in \textit{$<$interaction$>$.pot.in}. As already mentioned, \votca automatically creates potential tables to run a simulation. However, it does not know how to run a coarse-grained simulation. Therefore, all files needed to run a coarse-grained simulation, except for the potentials that are iteratively refined, must be provided and added to the \hyperlink{\cgref{inverse.filelist}}{filelist} in the settings \xml-file. If an atomistic topology and a mapping definition are present, \votca offers tools to assist the setup of a  coarse-grained topology (see chapter \ref{sec:usage:cgrun}).

To get an overview of how input files look like, it is suggested to take a look at one of the tutorials provided on \votcaweb.

In what follows we describe how to set up the iterative coarse-graining, run the main script, continue the run, and add customized scripts.

\subsection{Preparing the run}
\label{sec:preparing_the_run}
To start the first iteration, one has to prepare the input for the sampling program. This means that all files for running a coarse-grained simulation must be present and described in a separate \xml file, in our case \texttt{settings.xml} (see sec. \ref{sec:setting_files} for details). An extract from this file is given below. The only exception are tabulated potentials, which will be created and updated by the script in the course of the iterative process.

The input files include: target distributions, initial guess (optional) and a list of interactions to be iteratively refined. As a target distribution, any table file can be given (e.g. \gromacs output from \texttt{g\_rdf}). The program automatically takes care to resample the table to the correct grid spacing according to the options provided in \texttt{settings.xml}.

The initial guess is normally taken as a potential of mean force and is generated by Boltzmann-inversion of the corresponding distribution function. It is written in \texttt{step\_000/<name>.pot.new}. If you want to manually specify the initial guess for a specific interaction, write the potential table to a file called \texttt{<name>.pot.in} in the folder where you plan to run the iterative procedure.

A list of interactions to be iteratively refined has to be given in the options file. As an example, the \texttt{setting.xml} file for a propane is shown in listing~\ref{list:settings}. For more details,  see the full description of all options in ref.~\ref{sec:ref_options}.
\begin{figure}
\centering
\framebox[\textwidth]{\lstinputlisting{functionality/settings.xml}
}
\caption{\texttt{settings.xml} file specifies interactions to be refined, grid spacings, sampling engine, and the iterative method. The complete file can be found in the \texttt{propane/ibm} tutorial. 
\label{list:settings}
}
\end{figure}

\subsection{Starting the iterative process}
\label{sec:starting_iterative_process}
After all input files have been set up, the run can be started by
\begin{verbatim}
  csg_inverse --options settings.xml
\end{verbatim}

Each iteration is stored in a separate directory, named \texttt{step\_<iteration>}. \texttt{step\_000} is a special folder which contains the initial setup. For each new iteration, the files required to run the CG simulation (as specified in the config file) are copied to the current working directory. The updated potentials  are copied from the last step, \texttt{step\_<n-1>/<interaction>.pot.new}, and used as the new working potentials \texttt{step\_<n>/<interaction>.pot.cur}.

After the run preparation, all potentials are converted into the format of the sampling program and the simulation starts. Once the sampling has finished, analysis programs generate new distributions, which are stored in \texttt{<interaction>.dist.new}, and new potential updates, stored in \texttt{<interaction>.dpot.new}. 

Before adding the update to the old potential, it can be processed in the \texttt{post\_update} step. For each script that is specified in the postupdate, \texttt{<interaction>.dpot.new} is renamed  to \texttt{<interaction>.dpot.old} and stored in \texttt{<interaction>.dpot.<a-number>} before the processing script is called. Each processing script  uses the current potential update \texttt{<interaction>.dpot.cur} and writes the processed update to \texttt{<interaction>.dpot.new}. As an example, a pressure correction is implemented as a postupdate script within this framework.

After all postupdate scripts have been called, the update is added to the potential and the new potential \texttt{<interaction>.pot.new} is written. Additional post-processing of the potential can be performed in the \texttt{post\_add} step which is analogous to the \texttt{post\_update} step except for a potential instead of an update.

To summarize, we list all standard output files for each iterative step:

\begin{tabular}{ll}
\texttt{*.dist.new} & distribution functions of the current step \\
\texttt{*.dpot.new} & the final potential update, created by \texttt{calc\_update} \\
\texttt{*.dpot.<number>} & for each postupdate script, the \texttt{.dpot.new} is saved and a new one\\
&is created\\
\texttt{*.pot.cur} & the current potential used for the actual run\\
\texttt{*.pot.new} & the new potential after the add step \\
\texttt{*.pot.<number>} & same as \texttt{dpot.<number>} but for \texttt{post\_add}
\end{tabular}

If a sub-step fails during the iteration, additional information can be found in the log file. The name of the log file is specified in the steering \xml file.

\subsection{Restarting and continuing}
The interrupted or finished iterative process can be restarted either by extending a finished run or by restarting the interrupted run. When the script \prog{csg_inverse} is called, it automatically checks for a file called \texttt{done} in the current directory. If this file is found, the program assumes that the run is finished. To extend the run, simply increase \cgopt{inverse.iterations_max} in the settings file and remove the file called \texttt{done}. After that, \prog{csg_inverse} can be restarted, which will automatically recognize existing steps and continue after the last one.

If the iteration was interrupted, the script \prog{csg_inverse} might not be able to restart on its own. In this case, the easiest solution is to delete the last step and start again. The script will then repeat the last step and continue. However, this method is not always practical since sampling and analysis might be time-consuming and the run might have only crashed due to some inadequate post processing option. To avoid repeating the entire run, the script \prog{csg_inverse} creates a file with restart points and labels already completed steps such as simulation, analysis, etc. The file name is specified in the option \cgopt{inverse.restart_file}. If specific actions should be redone, one can simply remove the corresponding lines from this file. Note that a file \texttt{done} is also created in each folder for those steps which have been successfully finished.


% #####################################################################################################################
% #####################################################################################################################

\section{Iterative Boltzmann Inversion}
\subsection{Input preparation}
This section describes the usage of \ibi, implemented within the scripting framework described in the previous section \ref{sec:iterative_workflow}. It is suggested to get a basic understanding of this framework before proceeding.

\ibi so far only supports iterative refinement of non-bonded interactions. An outline of the workflow for performing \ibi is given in \fig{fig:flow_ibi}.

To specify Iterative Boltzmann Inversion as algorithm in the script, add \texttt{ibi} in the \texttt{method} section of the \xml setting file as shown below.

\begin{lstlisting}
  <cg>
    ...
    <inverse>
      <method>ibi</method>
    </inverse>
  </cg>
\end{lstlisting}


% #####################################################################################################################
% #####################################################################################################################

\section{Inverse Monte Carlo}
\label{sec:iterative_methods_imc}
In this section, additional options are described to run \imc coarse graining. The usage of \imc is similar to the one of \ibi and understanding the use of the scripting framework described in chapter~\ref{sec:iterative_workflow} is necessary.

\textbf{WARNING: multicomponent \imc is still experimental!}

\subsection{General considerations}
In comparison to \ibi, \imc needs significantly more statistics to calculate the potential update\cite{Ruehle:2009.a}. It is advisable to perform smoothing on the potential update. Smoothing can be performed as described in \sect{ref:ibi:optimize}. In addition, \imc can lead to problems related to finite size: for methanol, an undersized system proved to lead to a linear shift in the potential\cite{Ruehle:2009.a}. It is therefore always necessary to check that the system size is sufficiently large and that runlength csg smoothing iterations are well balanced.

\subsection{Additional mapping for statistics}
The program \prog{csg_stat} is used for evaluating the \imc matrix. Although the matrix only acts on the coarse-grained system here, it still needs a mapping file to work. This will improve with one of the next releases to simplify the setup. The mapping file needs to be a one to one mapping of the coarse grained system, e.g. for coarse graining \spce water, the mapping file looks as follows:
\begin{lstlisting}
  </cg_molecule>
    <name>SOL</name> 
    <ident>SOL</ident>
    <topology>
      <cg_beads>
        <cg_bead>
          <name>CG</name>
          <type>CG</type>
          <mapping>A</mapping>
          <beads>
            1:SOL:CG 
          </beads>
        </cg_bead>
      </cg_beads>
    </topology>
    <maps>
      <map>
        <name>A</name>
        <weights>1</weights>
      </map>
    </maps>
  </cg_molecule>
\end{lstlisting}

\subsection{Correlation groups}
Unlike \ibi, \imc also takes cross-correlations of interactions into account in order to calculate the update. However, it might not always be beneficial to evaluate cross-correlations of all pairs of interactions. By specifying \interopt{inverse.imc.group}, \votca allows to define groups of interactions, amongst which cross-correlations are taken into account, where \interopt{inverse.imc.group} can be any name.

\begin{lstlisting}
  <non-bonded>
    <name>CG-CG</name>
    <type1>CG</type1>
    <type2>CG</type2>
    ...
    <imc>
      <group>solvent</group>
   </imc>
  </non-bonded>
  <non-bonded>
\end{lstlisting}


% #####################################################################################################################
% #####################################################################################################################

\section{Pressure correction}

The pressure of the coarse-grained system usually does not match the pressure of the full atomistic system. This is because iterative Boltzmann inversion only targets structural properties but not thermodynamic properties. In order correct the pressure in such a way that it matches the target pressure (\interopt{inverse.p_target})., different strategies have been used based on small modifications of the potential. The correction can be enable by adding pressure to the list of \interopt{inverse.post_update} scripts. The type of pressure correction is selected by setting \interopt{inverse.post_update_options.pressure.type}.

\subsection{Simple pressure correction}
In ref.\cite{Reith:2003} a simple linear attractive potential was added to the coarse-grained potential
\begin{equation}
  \Delta V(r)=A \left( 1-\frac{r}{r_{cutoff}} \right) \,,
\end{equation}
with prefactor $A$
\begin{equation}
  A = -\sign(\Delta P)0.1k_{B}T\min(1,|f\Delta P) \,,
\end{equation}
$\Delta p=P_i-P_\text{target}$, and scaling factor $f$ and $P_\text{target}$ can be specified in the settings file as \interopt{inverse.post_update_options.pressure.simple.scale} and \interopt{inverse.p_target}.

As an example for a block doing simple pressure correction, every third interaction is
\begin{lstlisting}
<post_update>pressure</post_update>
<post_update_options>
  <pressure>
    <type>simple</type>
    <do>0 0 1</do>
    <simple>
      <scale>0.0003</scale>
    </simple>
  </pressure
</post_update_options>
\end{lstlisting}
Here, \interopt{inverse.post_update_options.pressure.simple.scale} is the scaling factor $f$. In order to get the correct pressure it can become necessary to tune the scaling factor $f$ during the iterative process.

\subsection{Advanced pressure correction}
In \cite{Wang:2009} a pressure correction based on the virial expression of the pressure was introduced. The potential term remains as in the simple form while a different sturcture of the $A$ factor is used:
\begin{equation}
  A = \left[\frac{-2\pi\rho^{2}}{3r_{cut}}\int_{0}^{r_{cut}}r^{3}g_{i}(r)dr\right]A_{i}=\Delta P.
\end{equation}
This factor requires the particle density $ \rho $ as additional input parameter, which is added as  \interopt{inverse.particle_dens} in the input file.

\section{Kirkwood-Buff correction}
In order to reproduce the exact Kirkwood-Buff ingetrals (KBIs), an correction term can be added into the coarse-grained potential~\cite{Ganguly:2012},
\begin{equation}
  \Delta U_{ij}^{(n)}(r) = \frac{k_{B}T}\;A\;(G_{ij}^{(n)} - G_{ij}^\text{ref})\left(1- \frac{r}{r_\text{ramp}}\right),
\end{equation}
where $G_{ij}^{(ref)}$ is the KBI calculated from the reference all-atom simulation and $G_{ij}^{(n)}$ is the KBI 
after the $n^{th}$ iteration.

The Kirkwood-Buff integrals are calculated from the radial distribution functions as follows:
\begin{equation}
G_{ij} = 4\pi \int_0^\infty \left[ g_{ij}(r) - 1\right] r^2 dr~.
\label{eq:kbi}
\end{equation}
For simulations of finite box size we calculate the running integral up to distance $R$
\begin{equation}
  G_{ij}(R) = 4\pi \int_0^R \left[ g_{ij}(r) - 1\right] r^2 dr~.
\end{equation}
The average of those running integrals in the interval, where $G_{ij}(R)$ gets flat, gives a good estimate for $G_{ij}$:
\begin{equation}
  G_{ij}\approx<G_{ij}(R)>|_{R=r_1}^{R=r_2}
\end{equation}
As an example for a block doing Kirkwood-Buff correction, every iteraction without doing potential update
\begin{lstlisting}
<do_potential>0</do_potential>
<post_update>kbibi</post_update>
<post_update_options>
  <kbibi>
    <do>1</do>
    <start>1.0</start>
    <stop>1.4</stop>
    <factor>0.05</factor>
    <r_ramp>1.4</r_ramp>
  </kbibi>
</post_update_options>
\end{lstlisting}
Here, \interopt{inverse.post_update_options.kbibi.factor} is the scaling factor $A$. \interopt{inverse.post_update_options.kbibi.start} is $r_1$ and \interopt{inverse.post_update_options.kbibi.stop} is $r_2$ used to calculate the average of $G_{ij}(R)$.
\section{Runtime optimization}
\label{ref:ibi:optimize}
Most time per iteration is spent on running the coarse-grained system and on calculating the statistics. To get a feeling on how much statistics is needed, it is recommended to plot the distribution functions and check whether they are sufficiently smooth. Bad statistics lead to rough potential updates which might cause the iterative refinement to fail. All runs should be long enough to produce distributions/rdfs of reasonable quality.

Often, runtime can be improved by smoothing the potential updates. Our experience has shown that it is better to smooth the potential update instead of the rdf or potential itself. If the potential or rdf is smoothed, sharp features like the first peak in \spce water might get lost. Smoothing on the delta potential works quite well, since the sharp features are already present from the initial guess. By applying iterations of a simple triangular smoothing ($ \Delta U_i = 0.25 \Delta U_{i-1} + 0.5\Delta U_i + 0.25\Delta U_{i+1} $), a reasonable coarse-grained potential for \spce water could be produced in less than 10 minutes. Smoothing is implemented as a post\_update script and can be enabled by adding
\begin{lstlisting}
  <post_update>smooth</post_update>
  <post_update_options>
    <smooth>
        <iterations>2</iterations>
    </smooth>
  </post_update_options>
\end{lstlisting}
to the inverse section of an interaction in the settings \xml file.


\section{Thermodynamic force}
The thermodynamic force method is an iterative procedure to determine an external field that can correct for density variations. This has been prooven to be usefull for multi-scale simulations where all-atom and coarse-grained representations are simulated concurrently in one simulation. The AdresS simulation scheme provides a protocol for such simulations.\\
The thermodynamic force is updated from the density profile in each simulation step as:
\begin{equation}
  {\bf F}_\text{th}^{i+1}({\bf r}) = {\bf F}_\text{th}^i({\bf r}) - \frac{m_\text{mol}}{\rho_0^2 \kappa_T^\text{at}}\nabla \rho_i({\bf r})
\label{iterative_force_grad_dens}
\end{equation}
where $\rho_i({\bf r})$ can be either a density along one of the box axis or a radial density, this is specified by the adress paramater adress\_type in the gromacs mdp file.
In order to use the thermodynamic force iteration, VOTCA must be used together with the gromacs-4.6.1 or later. To check whether your gromacs version support this type \begin{verbatim}mdrun -h \end{verbatim} and look for the -tabletf option.
A tutorial simulation set can be found in the tutorials (spce/tf) which performs the thermodynamic force iteration for spc/e water coupled to a coarse-grained spc/e water.\\
The method is selected by specifying
\begin{lstlisting}
<method>tf</method>
\end{lstlisting}
in the inverse section.
For each interaction type additional options have to be specified in the settings.xml file. To specify in which region the thermodynamic force should be nonzero, the min and max properties are used.
A smoothing function proportional to $cos^2(r)$ is used to make the force go smoothly to zero at the region specified by min and max.
Additonally a 'tf' section is needed for each interaction type
\begin{lstlisting}
<non-bonded>
    <name>SOL</name>
    <min>1.4</min>
    <max>3.1</max>
   <step>0.01</step>
    <tf>
	<spline_start>0.9</spline_start>
	<spline_end>3.6</spline_end>
	<spline_step>0.4</spline_step>
	<molname>SOL</molname>
	<prefactor>0.01382</prefactor>
    </tf>
    <inverse>
	<target>dens.SOL.xvg</target>
	(...)
    </inverse>
</non-bonded>
\end{lstlisting}
Usually the density profile fluctuates too much to obtain a force directly from the gradient. Thus spline interpolation is used to smooth the force. To specify the spline interpolation range the spline\_start and spline\_end parameters are used. These can define a larger region than between min and max as it is sometimes usefull to extend the spline fit for numerical stability. The parameter spline\_step sepcifies the bin width of the fit grid (see \prog{csg_resample} for more). 
The field 'molname' specifies the molecule (as defined in the gromacs topology) used for calculating the density. 
The prefactor $\frac{m_\text{mol}}{\rho_0^2 \kappa_T^\text{at}}$ appearing in eq \ref{iterative_force_grad_dens} is specified in the 'prefactor' field. Its units are $\text{kJ\;nm}^3/(\text{mol\;u})$ in the case of \gromacs due to the fact $\rho$ is calculated in $\text{u}/\text{nm}^3$ and force units are $\text{kJ}/(\text{mol\;nm})$. If two representations of the fluid with different compressibilities are coupled, it is advisable to use the higher compressibility in order to stabilize the iteration.
A target density file has to be specified for each interaction type, in most cases this will containt a flat density profile at the equilibrium density $\rho_0$. \\
